# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
# install.packages(c("tidyverse", "estimatr", "car"))

library(tidyverse)
library(estimatr)
# car called with car::recode() but not loaded because of conflicts with dplyr::recode()

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

dv_df_mturk <- dv_df
dv_df_elite <- filter(dv_df, dv_cat != "climate")

all_fits_mturk <- NULL
for(i in 1:nrow(dv_df_mturk)) {

    fit_df_w1 <-
    lm_robust(formula(paste0(dv_df_mturk$dvs[i], "_w1 ~ Z")),
              data = mturk_opeds) %>%
    tidy %>%
    mutate(wave = "wave_1", dvs = dv_df_mturk$dvs[i])
  
  fit_df_w2 <-
    lm_robust(formula(paste0(dv_df_mturk$dvs[i], "_w2 ~ Z")),
              data = mturk_opeds) %>%
    tidy %>%
    mutate(wave = "wave_2", dvs = dv_df_mturk$dvs[i])
  
  fit_df_w3 <-
    lm_robust(formula(paste0(dv_df_mturk$dvs[i], "_w3 ~ Z")),
              data = mturk_opeds) %>%
    tidy %>%
    mutate(wave = "wave_3", dvs = dv_df_mturk$dvs[i])
  
  all_fits_mturk <-
    bind_rows(all_fits_mturk, fit_df_w1, fit_df_w2, fit_df_w3)
}

all_fits_elite <- NULL
for (i in 1:nrow(dv_df_elite)) {
  fit_df_w1 <-
    lm_robust(formula(paste0(dv_df_elite$dvs[i], "_w1 ~ Z")),
              data = elite_opeds) %>%
    tidy %>%
    mutate(wave = "wave_1", dvs = dv_df_elite$dvs[i])
  
  fit_df_w2 <-
    lm_robust(formula(paste0(dv_df_elite$dvs[i], "_w2 ~ Z")),
              data = elite_opeds) %>%
    tidy %>%
    mutate(wave = "wave_2", dvs = dv_df_elite$dvs[i])
  
  all_fits_elite <- bind_rows(all_fits_elite, fit_df_w1, fit_df_w2)
}

all_fits <- bind_rows(elite = all_fits_elite, mturk = all_fits_mturk, .id = "study")

all_fits <-
  all_fits %>%
  filter(term != "(Intercept)") %>%
  mutate(stars = car::recode(p.value, "0:0.001 = '***';0.001:0.01 = '**'; 0.01:0.05 = '*'; 0.05:1 = ''"),
         treatment = car::recode(term, "'Zamtrak' = 'Amtrak'; 'Zclimate'='Climate'; 'Zpaul'='Paul'; 'Zveterans'='Veterans'; 'Zwallstreet'='Wall Street'"),
         treatment = factor(treatment, levels = c("Climate", "Wall Street", "Amtrak", "Veterans", "Paul"))) %>%
  left_join(dv_df)


g1 <-
  filter(all_fits, study == "mturk", wave == "wave_1") %>%
  ggplot(aes(x = treatment, y = dv_labels)) +
  geom_tile(aes(fill = estimate)) +
  geom_text(aes(label = stars)) +
  scale_y_discrete(drop = FALSE) +
  scale_fill_gradient2(
    low = "#205C8A",
    mid = "white",
    high = "#C67800",
    midpoint = 0,
    limits = c(-1, 1)
  ) +
  xlab("Op-Ed Shown") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    legend.key.height = unit(4, "lines")
  )


g2 <- g1 %+% filter(all_fits, study == "mturk", wave == "wave_2")
g3 <- g1 %+% filter(all_fits, study == "mturk", wave == "wave_3")
g4 <- g1 %+% filter(all_fits, study == "elite", wave == "wave_1")
g5 <- g1 %+% filter(all_fits, study == "elite", wave == "wave_2")

g1
g2
g3
g4
g5




